---
title: "Camelid size"
author: "Juliana Rubinatto Serrano"
date: "2024-01-29"
output: html_document
---

An R vignette for analyzing South American camelid size and taxonomy based on phalanx measurement. This work is based on the R research compendium used in the statistical analysis and visualization for:

deFrance, S. D y Juliana Rubinatto Serrano (2025) "Camélidos Tihuanacos en la Colonia de Moquegua al Sur del Perú: Variaciones de Tamaño y Morfotipos". *Chungara: Humanos y camélidos: Interacciones sociales e historia evolutiva.*

# Data and set up packages

```{r setup}
library(tidyr)
library(dplyr)
library(reshape2)
library(zoolog)
library(ggplot2)
library(ggpubr)
library(rstatix)
```

```{r data wrangling}
data <- read.csv("data/data.csv")

```

#Exploratory Data Analysis and initial graphics

```{r EDA graphics}

data$species <- factor(data$species,levels = c("Moquegua camélido (n = 298)", "Le Neün et al. (2023) camélido arqueologico (n = 125)", "Alpaca moderna (n = 22)", "Guanaco moderno (n = 179)", "Llama moderna (n = 109)", "Vicuña moderna (n = 109)"))

# Displaying all data from all samples
## used to generate Figure 3 in deFrance and Rubinatto Serrano (2025)
plot_by_sample <- ggplot(data, aes(x=H, y=Bp, 
                                 shape = species, 
                                 color = species)) + 
  geom_point(size = 10) + 
  annotate(geom = "text",x=13, y=10, label="Pequeño", size = 15) +
  annotate(geom = "text",x=25, y=10, label="Grande", size = 15) +
  scale_color_manual(values=c("#7851A9", "#FF7900", "#DF30FF",
                              "#48ADDE", "#A5DE00", "#FF9111"),
                     labels = c("Moquegua camélido (n = 298)", 
                                "Le Neün et al. (2023) \ncamélido arqueologico (n = 125)",
                                "Alpaca moderna (n = 22)", 
                                "Guanaco moderno (n = 179)", 
                                "Llama moderna (n = 109)", 
                                "Vicuña moderna (n = 109)")) +
  scale_shape_manual(values=c(19, 20, 18, 18, 18, 18),
                     labels = c("Moquegua camélido (n = 298)", 
                                "Le Neün et al. (2023) \ncamélido arqueologico (n = 125)",
                                "Alpaca moderna (n = 22)", 
                                "Guanaco moderno (n = 179)", 
                                "Llama moderna (n = 109)", 
                                "Vicuña moderna (n = 109)"))
  
plot_by_sample + labs(x = "V3 (mm)", y = "V2 (mm)", color = "", shape = "") + 
  theme(legend.position = "bottom",
        axis.title = element_text(size = 35), 
        plot.title = element_text(size = 35), 
        legend.text = element_text(size = 40))+
  geom_abline(intercept = 45.5, slope = -1.5, color="red", 
                 linetype="dotdash", size=.5) +
  geom_abline(intercept = 42, slope = -1.5, color="red", 
                 linetype="twodash", size=.5)


# Displaying data from one of the groups by archaeological site
## used to generate Figure 2 in deFrance and Rubinatto Serrano (2025)
Moquegua <- data %>% filter(group == "Moquegua (n = 298)")
long_data_Moquegua <- Moquegua %>% select(Sitio, Bp, H) %>% 
  gather(key = "variable", value = "value", -Sitio) #long format
long_data_Moquegua_GL <- Moquegua %>% select(Sitio, GL) %>% 
  gather(key = "variable", value = "value", -Sitio)

species <- data %>% filter(group != "Le Neün et al. (2023) arqueologica (n = 125)")

species$species <- factor(species$species,levels = c("camélido", "alpaca", 
                                                     "vicuña", "llama", "guanaco"))
long_data_species <- species %>% select(species, Bp, H) %>% 
  gather(key = "variable", value = "value", -species) #long format
long_data_speciesGL <- species %>% select(species, GL) %>% 
  gather(key = "variable", value = "value", -species) #long format

##Scatterplot by site

plot_by_site <- ggplot(Moquegua, aes(x=H, y=Bp, color= Sitio)) + 
  geom_point(size = 10) + 
  annotate(geom = "text",x=13, y=10, label="Pequeño", size = 15) +
  annotate(geom = "text",x=25, y=10, label="Grande", size = 15) +
  scale_color_manual(values=c("#DF30FF","#48ADDE", "#A5DE00", "#FF9111"),
                     labels=c("M1 - Chen Chen\n (n = 22)", "M10 - Omo\n (n = 71)",
                              "M12 - Omo\n (n = 8)","M43 - Rio Muerto\n (n = 197)"))

plot_by_site + labs(x = "V3 (mm)", y = "V2 (mm)") +
  theme(legend.position = "bottom", axis.title = element_text(size = 35),
        legend.title = element_blank(),
        plot.title = element_text(size = 35), legend.text = element_text(size = 40)) +
  geom_abline(intercept = 45.5, slope = -1.5, color="red", 
                 linetype="dotdash", size=.5) +
  geom_abline(intercept = 42, slope = -1.5, color="red", 
                 linetype="twodash", size=.5)


#Boxplot sites and samples
## used to generate Figure 4 in deFrance and Rubinatto Serrano (2025)

###By site
boxplot1 <- ggplot(long_data_Moquegua, aes(x = variable, y=value, color = Sitio)) + 
  geom_boxplot(linewidth = 2) +
   scale_color_manual(values=c("#DF30FF","#48ADDE", "#A5DE00", "#FF9111"),
                     labels=c("M1 - Chen Chen\n (n = 22)", "M10 - Omo\n (n = 71)",
                              "M12 - Omo\n (n = 8)","M43 - Rio Muerto\n (n = 197)")) + 
  labs(x = "Bp = V2 y H = V3", y = "Valores (mm)") +
  theme(legend.position = "none", legend.title = element_blank(),
        axis.title = element_text(size = 35),
        plot.title = element_text(size = 35), legend.text = element_text(size = 40))

boxplot2 <- ggplot(long_data_Moquegua_GL, aes(x = variable, y=value, color = Sitio)) + 
  geom_boxplot(linewidth = 2) +
   scale_color_manual(values=c("#DF30FF","#48ADDE", "#A5DE00", "#FF9111"),
                     labels=c("M1 - Chen Chen\n (n = 22)", "M10 - Omo\n (n = 71)",
                              "M12 - Omo\n (n = 8)","M43 - Rio Muerto\n (n = 197)")) + 
  labs(x = "GL = V1", y = "") +
  theme(legend.position = "none", axis.title = element_text(size = 35))

###By sample
long_data <- data %>% select(group, Bp, H) %>% 
  gather(key = "variable", value = "value", -group) #long format
long_dataGL <- data %>% select(group, GL) %>% 
  gather(key = "variable", value = "value", -group)

boxplot3 <- ggplot(long_data, aes(x = variable, y=value, color = group)) + 
  geom_boxplot(linewidth = 2) +
   scale_color_manual(values=c("#43984A", "#FF7900", "#7851A9"),
                      labels=c("Le Neün et al. (2023)\narqueologica (n = 125)", 
                              "Le Neün et al. (2023)\nmoderna (n = 419)",
                              "Moquegua \n(n = 298)")) + 
  labs(x = "Bp = V2 y H = V3", y = "Valores (mm)") +
  theme(legend.position = "none", legend.title = element_blank(),
        axis.title = element_text(size = 35), 
        plot.title = element_text(size = 35), legend.text = element_text(size = 40))

boxplot4 <- ggplot(long_dataGL, aes(x = variable, y=value, color = group)) + 
  geom_boxplot(linewidth = 2) +
   scale_color_manual(values=c("#43984A", "#FF7900", "#7851A9"),
                      labels=c("Le Neün et al. (2023)\narqueologica (n = 125)", 
                              "Le Neün et al. (2023)\nmoderna (n = 419)",
                              "Moquegua \n(n = 298)")) + 
  labs(x = "GL = V1", y = "") +
  theme(legend.position = "none", axis.title = element_text(size = 35))



###Combining boxplots
big1 <- ggarrange(boxplot1, boxplot2, common.legend = TRUE, legend = "bottom")
big2 <- ggarrange(boxplot3, boxplot4, common.legend = TRUE, legend = "bottom")

ggarrange(big1, big2, nrow = 2, legend = "bottom", labels ="AUTO", font.label = list(size= 40))

```

# Statistical tests (ANOVA, MANOVA, and t-tests)

```{r statistical tests}
#Between samples
manova_samples <- manova(cbind(GL,Bp, H) ~ group, data = data) 
summary(manova_samples)
manova_samples1 <- manova(cbind(Bp, H) ~ group, data = data) 
summary(manova_samples1)

  ##geometric mean
  data_noNA <- subset(data, !is.na(mean_geo))
  anova_geo <- oneway.test(mean_geo ~ group, data = data_noNA, var.equal = TRUE)
  anova_geo
  
  GL <- data %>% select(group, GL)
  GL_noNA <- subset(data, !is.na(GL))
  res.aov <- oneway.test(GL ~ group, data = GL_noNA, var.equal = TRUE)
res.aov

  ##two-sided t-tests to compare

    ###Moquegua y le Neun et al (2023) arqueologica
    archaeological <- data %>% filter(group == "Moquegua (n = 298)" | 
                                        group == "Le Neün et al. (2023) arqueologica (n = 125)")
    arqueo <- t.test(cbind(Bp, H) ~ factor(archaeological$group), data = archaeological, 
                     alternative = "less", var.equal = TRUE)
    
    ###Moquegua y le Neun et al (2023) moderna
    mixed <- data %>% filter(group == "Moquegua (n = 298)" | 
                                        group == "Le Neün et al. (2023) moderna (n = 419)")
    mixed_test <- t.test(cbind(Bp, H) ~ factor(mixed$group), data = mixed, 
                     alternative = "greater", var.equal = TRUE)
     
    ###le Neun et al (2023) moderna y arqueologica
    leNeunetal <- data %>% filter(group == "Le Neün et al. (2023) arqueologica (n = 125)" | 
                                        group == "Le Neün et al. (2023) moderna (n = 419)")
    leNeunetal_test <- t.test(cbind(Bp, H) ~ factor(leNeunetal$group), data = leNeunetal, 
                     alternative = "less", var.equal = TRUE)
    
    print(arqueo)
    print(mixed_test)
    print(leNeunetal_test)


#Tiwanaku sites in Moquegua
manova_sites <- manova(cbind(GL, Bp, H) ~ Sitio, data = Moquegua) 
summary(manova_sites)
    
manova_sites1 <- manova(cbind(Bp, H) ~ Sitio, data = Moquegua) 
summary(manova_sites1)


#Fore and Hindlimb
  ##forelimb
  just_fore <- data %>% filter(fore_hind == "fore")
  manova_fore <- manova(cbind(GL, Bp, H) ~ group, data = just_fore) 
  summary(manova_fore) 
  
  ##hindlimb
  just_hind <- data %>% filter(fore_hind == "hind")
  manova_hind <- manova(cbind(GL, Bp, H) ~ group, data = just_hind) 
  summary(manova_hind) 


#Between samples of Peru and Bolivia
  Bolivia <- data %>% 
    filter(Origen == "Khonkho Wankane " | Origen == "Chiripa" | Origen == "Bolivia") %>%
    mutate(country = "Bolivia")
  Peru <- data %>% filter(source != "Moquegua" & Origen == "Peru") %>% mutate(country = "Peru")
  Moquegua2 <- Moquegua %>% mutate(country = "Tiwanaku")
  
  Peru_Bolivia <- bind_rows(Bolivia, Peru, Moquegua2)
  
  manova_origin <- manova(cbind(GL, Bp, H) ~ country, data = Peru_Bolivia) 
  summary(manova_origin) 

    ##Peru samples
    Peru_t <- t.test(cbind(Bp, H) ~ country, 
                     data = (filter(Peru_Bolivia, country == "Peru" | country == "Tiwanaku")), 
                     alternative = "less", var.equal = FALSE)
    print(Peru_t)
    
    ##Moquegua and Bolivia
    Bolivia_Tiwanaku <- t.test(cbind(GL, Bp, H) ~ country, 
                               data = (filter(Peru_Bolivia, 
                                              country == "Bolivia" | country == "Tiwanaku")), 
                               alternative = "greater", var.equal = FALSE)
    print(Bolivia_Tiwanaku)
    
    ##Bolivia and Peru
    Bolivia_Peru <- t.test(cbind(Bp, H) ~ country, 
                           data = (filter(Peru_Bolivia, country == "Bolivia" | country == "Peru")), 
                           alternative = "less", var.equal = FALSE)
    print(Bolivia_Peru)
```

# Visualization of comparison of samples from Bolivia and Peru

```{r visualization regional origin}

#Boxplot regional origin
## used to generate Figure 5 in deFrance and Rubinatto Serrano (2025)

Peru_Bolivia$country <- factor(Peru_Bolivia$country, 
                               levels = c("Peru", "Bolivia", "Tiwanaku"))
  Peru_Bolivia_long <- Peru_Bolivia %>% 
    select(Bp, H, country) %>% 
    gather(key = "variable", value = "value", -country)
  Peru_Bolivia_long_GL <- Peru_Bolivia %>% 
    select(GL, country) %>% 
    gather(key = "variable", value = "value", -country)
  
  boxplot5 <- ggplot(Peru_Bolivia_long, aes(x = variable, y=value, color = country)) + 
    geom_boxplot(linewidth = 3) +
    scale_color_manual(values=c("#FDB147", "#11574A" , "#7851A9"), 
                       labels=c("Le Neün et al. (2023)\nPeru (n = 136)", 
                                "Le Neün et al. (2023)\nBolivia (n = 46)",
                                "Moquegua (n = 298)")) + 
    labs(x = "Bp = V2 y H = V3", y = "Valores (mm)") +
    theme(legend.position = "bottom", legend.title = element_blank(),
          axis.title = element_text(size = 35), 
          plot.title = element_text(size = 35), legend.text = element_text(size = 40))
  
  boxplot6 <- ggplot(Peru_Bolivia_long_GL, aes(x = variable, y=value, color = country)) + 
  geom_boxplot(linewidth = 3) +
   scale_color_manual(values=c("#FDB147", "#11574A" , "#7851A9"), 
                       labels=c("Le Neün et al. (2023)\nPeru (n = 136)", 
                                "Le Neün et al. (2023)\nBolivia (n = 46)",
                                "Moquegua (n = 298)")) + 
    labs(x = "GL = V1", y = "") +
    theme(legend.position = "none", axis.title = element_text(size = 35))
  
ggarrange(boxplot5, boxplot6, common.legend = TRUE, legend = "bottom")
  
#Boxplot Bolivia and Moquegua
## used to generate Figure 6 in deFrance and Rubinatto Serrano (2025)
  Bolivia_Tiwanaku_long <- Peru_Bolivia %>% 
    filter(country == "Tiwanaku" | country == "Bolivia") %>% 
    select(Bp, H, comparison) %>% 
    gather(key = "variable", value = "value", -comparison)
  Bolivia_Tiwanaku_long_GL <- Peru_Bolivia %>% 
    filter(country == "Tiwanaku" | country == "Bolivia") %>% 
    select(GL, comparison) %>% 
    gather(key = "variable", value = "value", -comparison)
  
  boxplot7 <- ggplot(Bolivia_Tiwanaku_long, aes(x = variable, y=value, color = comparison)) + 
    geom_boxplot(linewidth = 3) +
    scale_color_manual(values=c("#FDB147", "#ADB147", "#71874A" , "#7851A9")) +
    labs(x = "Bp = V2 y H = V3", y = "Valores (mm)") +
    theme(legend.position = "bottom", legend.title = element_blank(),
          axis.title = element_text(size = 35), 
          plot.title = element_text(size = 35), legend.text = element_text(size = 40))
  
  boxplot8 <- ggplot(Bolivia_Tiwanaku_long_GL, aes(x = variable, y=value, color = comparison)) + 
  geom_boxplot(linewidth = 3) +
   scale_color_manual(values=c("#FDB147", "#ADB147", "#71874A" , "#7851A9")) + 
    labs(x = "GL = V1", y = "") +
    theme(legend.position = "none", axis.title = element_text(size = 35))
  
ggarrange(boxplot7, boxplot8, common.legend = TRUE, legend = "bottom")
```

# Log size index (LSI) using zoolog package

```{r setting up for analysis}
Measurement <- data.frame(data)

#some changes to the original data frame are required to use the zoolog thesaurus 
Measurement <- select(Measurement, -14)
Measurement$EL <- "PH1" #creating the element 1st phalanx 
Measurement$Species <- "Camelid" #creating the species name Camelid

#because we will use a different reference specimen for the fore and hind limb phalanx, we will use separate dataframes for each
Moqu_fore <- filter(Measurement, group == "Moquegua (n = 298)" & fore_hind == "fore")
Moqu_hind <- filter(Measurement, group == "Moquegua (n = 298)" & fore_hind == "hind")

#THESAURUS UPDATE NOTE: because the original zoolog package does not include Camelid as a unique taxon, Cervus was arbitrarily chosen to add Camelid into the thesaurus. Since we are using custom reference values, having Camelid under Cervis will not affect the results but allow us to use the package.
taxonThesaurus <- AddToThesaurus(
  thesaurus = zoologThesaurus$taxon,
  newName = c("Camelid"), 
  category = "Cervus")
zoologThesaurus$taxon <- taxonThesaurus
WriteThesaurusSet(thesaurus = zoologThesaurus, file = "myThesaurus")

#Creating the reference average fore-hind and adding it to zoolog
reference1_Bp <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("Bp"),
                  Standard = c("18.89"))
reference1_Bp$Standard <- as.numeric(reference1_Bp$Standard)
reference1_H <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("H"),
                  Standard = c("17.96"))
reference1_H$Standard <- as.numeric(reference1_H$Standard)
reference1_Gl <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("GL"),
                  Standard = c("61.67"))
reference1_Gl$Standard <- as.numeric(reference1_Gl$Standard)
reference1_camelid = list(reference1_Bp, reference1_H, reference1_Gl)

#Creating the reference for forelimb phalanx and adding it to zoolog
reference2_Bp <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("Bp"),
                  Standard = c("19.92"))
reference2_Bp$Standard <- as.numeric(reference2_Bp$Standard)
reference2_H <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("H"),
                  Standard = c("18.6"))
reference2_H$Standard <- as.numeric(reference2_H$Standard)
reference2_Gl <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("GL"),
                  Standard = c("65.67"))
reference2_Gl$Standard <- as.numeric(reference2_Gl$Standard)
reference2_camelid = list(reference2_Bp, reference2_H, reference2_Gl)

#Creating the reference for forelimb phalanx and adding it to zoolog
reference3_Bp <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("Bp"),
                  Standard = c("17.86"))
reference3_Bp$Standard <- as.numeric(reference3_Bp$Standard)
reference3_H <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("H"),
                  Standard = c("17.32"))
reference3_H$Standard <- as.numeric(reference3_H$Standard)
reference3_Gl <- data.frame (TAX  = c("Cervus"),
                  EL = c("PH1"),
                  Measure = c("GL"),
                  Standard = c("57.68"))
reference3_Gl$Standard <- as.numeric(reference3_Gl$Standard)
reference3_camelid = list(reference3_Bp, reference3_H, reference3_Gl)

```

## Running reference 1 for all samples

```{r reference 1, warning = FALSE}
#calculate LSI for each variable
testLogsBp <- LogRatios(Measurement, ref = reference1_Bp)
testLogsH  <- LogRatios(Measurement, ref = reference1_H)
testLogsGl <- LogRatios(Measurement, ref = reference1_Gl)
#add these LSI values to the data frame
Logs <- cbind(testLogsBp, testLogsGl$logGL, testLogsH$logH)
colnames(Logs)[17] <- "logGL"
colnames(Logs)[18] <- "logH"
#Long format of relevant data for visualization
Logs_long <- Logs %>% 
    select(group, logBp, logGL, logH) %>% 
    gather(key = "variable", value = "value", -group)

#Boxplot
## used to generate Figure 7 in deFrance and Rubinatto Serrano (2025)
boxplot_logs <- ggplot(Logs_long, aes(x = variable, y = value, color = group)) + 
  geom_boxplot(linewidth = 2) + coord_flip() + ylim(-0.3, 0.2) +
  scale_color_manual(values=c("#43984A", "#FF7900", "#7851A9"),
                     labels=c("Le Neün et al. (2023)\narqueologica (n = 125)", 
                              "Le Neün et al. (2023)\nmoderna (n = 419)",
                              "Moquegua \n(n = 298 por Bp y H)\n(n = 251 por GL)")) 
 
boxplot_logs + labs(x = "Diferentes medidas", y = "Relación logarítmica", color = "") +
  geom_hline(yintercept = 0, colour="red", linetype = "dashed", linewidth = 1) + #line for reference specimen value
  theme(legend.position = "bottom", legend.title = element_blank(),
        axis.title = element_text(size = 35), plot.title = element_text(size = 35), 
        legend.text = element_text(size = 40))
```

## Running reference 2 and 3 for forelimbs and hindlimbs respectively

```{r reference 2 and 3}

#Forelimb
testLogsBp2 <- LogRatios(Moqu_fore, ref = reference2_Bp)
testLogsH2 <- LogRatios(Moqu_fore, ref = reference2_H)
testLogsGl2 <- LogRatios(Moqu_fore, ref = reference2_Gl)
Logs_fore <- cbind(testLogsBp2, testLogsGl2$logGL, testLogsH2$logH)
colnames(Logs_fore)[17] <- "logGL"
colnames(Logs_fore)[18] <- "logH"

#Hindlimb
testLogsBp3 <- LogRatios(Moqu_hind, ref = reference3_Bp)
testLogsH3 <- LogRatios(Moqu_hind, ref = reference3_H)
testLogsGl3 <- LogRatios(Moqu_hind, ref = reference3_Gl)
Logs_hind <- cbind(testLogsBp3, testLogsGl3$logGL, testLogsH3$logH)
colnames(Logs_hind)[17] <- "logGL"
colnames(Logs_hind)[18] <- "logH"

#Combining results
Logs_all <- bind_rows(Logs_fore, Logs_hind)
Logs_all_long <- Logs_all %>% 
  select(fore_hind, logBp, logGL, logH) %>% 
  gather(key = "variable", value = "value", -fore_hind)

#used to generate Figure 8 in deFrance and Rubinatto Serrano (2025)
supps.labs <- c("Falange delantera (n = 6)", "Falange trasera (n = 8)")
names(supps.labs) <- c("fore", "hind")

ggplot(Logs_all_long, aes(x = variable, y = value, 
                          color = fore_hind, shape = fore_hind)) +
  geom_jitter(width = 0, height = 0, size = 5) +
  facet_grid(~ fore_hind, 
             labeller = labeller(fore_hind = supps.labs)) +
  theme(legend.position="none", 
        axis.title = element_text(size = 35), 
        plot.title = element_text(size = 35)) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = 2, color = "red", linewidth = 1) + 
  scale_color_manual(values = c("#5851A9", "#01954A"), 
                     labels = c("Falange delantera", "Falange trasera")) +
  scale_shape_manual(values = c(17, 16), 
                     labels = c("Falange delantera", "Falange trasera")) +
  labs(x = "Relaciónes logarítmicas", y = "Valores") 

```


## Other elements

```{r}
others <- read_excel("data/others.xlsx")
others_long <- others %>% 
  select(Element, LSI) %>%
  gather(key = "variable", value = "value", -Element)

#used to generate Figure 9 in deFrance and Rubinatto Serrano (2025)
ggplot(others_long, aes(x = variable, y = value, color = Element, shape = Element)) +
  geom_jitter(width = 0, height = 0, size = 3) +
  facet_grid(~ Element) +
  theme(legend.position="none", axis.title = element_text(size = 14),
        axis.text.x = element_text(angle = 90)) +
  coord_flip() + scale_y_continuous(breaks=seq(0, 0.08, 0.04)) +
  geom_hline(yintercept = 0, linetype = 2, color = "red", linewidth = 1) + 
  labs(x = "Relaciónes logarítmicas", y = "Values") 
```
